Supongamos que queremos minimizar una función $f:\mathbb{R}\rightarrow \mathbb{R}$ y tenemos un kernel de transición $Q$ para realizar una búsqueda aleatoria de este mínimo.
1) Sea $x_0$ un candidato inicial del mínimo de $f$ y $y_0=f(x_0)$
2) Para $i$ desde $1$ hasta $N$
2.1) Simule $x_n\sim Q(\cdot|x_{n-1})$ y tome $y_n=f(x_n)$
3) Tome $y^\star=\min\left\{y_0,\ldots , y_N\right\}$ como aproximación al valor que minimiza $f$.
El algoritmo de recocimiento simulado genera una cadena en la cual se acepta la propuesta cuando $f(x_{\text{prop}}) \leq f(x_{\text{prev}}) $ En caso de que $f(x_{\text{prop}}) > f(x_{\text{prev}})$ se considera un paso de aceptación tipo Metropolis-Hastings haciendo uso de la función de densidad de probabilidad
$$ p(x)=\frac{e^{-\frac{f(x)}{T}}}{c_p} $$donde $c_p$ es la constante de normalización y el parámetro $T$ tiene la interpretación de temperatura en el contexto de física estadística; y una probabilidad de aceptación $$ \alpha(x_{\text{prop}}|x_{\text{prev}})=\frac{ p(x_{\text{prop}}) }{ p(x_{\text{prev}}) }= e^{-\frac{\left( f(x_{\text{prop}})-f(x_{\text{prev}})\right)}{T}} \mathbb{1}_{\left\{ f(x_{\text{prop}}) > f(x_{\text{prev}}) \right\}}. $$
Por lo que con probabilidad positiva en el caso de que $ f(x_{\text{prop}}) > f(x_{\text{prev}}$ podemos aceptar $x_{\text{prop}}$ a pesar de que este valor es un peor candidato para la minimización de $f$. Observemos que si la temperatura $T$ es cercana a cero entnces la probabilidad de aceptación $\alpha$ será cercana a cero. Esto puede ser de utilidad para salir de regiones cercanas a un mínimo local cuando la temperatura $T$ no es cercana a cero.
1) Sea $x_0$ un candidato inicial del mínimo de $f$ y $y_0=f(x_0)$
2) Para $i$ desde $1$ hasta $N$
2.1) Simule $x_{\text{prop}}\sim Q(\cdot|x_{n-1})$ y tome $y_{\text{prop}}=f(x_n)$
2.2) Si $y_{\text{prop}} \leq y_{i-1}$ tome $x_i=x_{\text{prop}}$ y $y_i=y_{\text{prop}}$. Proceda a 2.3) sólo en el caso contrario $y_{\text{prop}} > y_{i-1}$
2.3) Tome $U\sim \text{Uniforme}(0,1)$
2.3.1) Si $U < \alpha(x_{\text{prop}}|x_{i-1})$ entonces tome $x_i=x_{\text{prop}}$ y $y_i=y_{\text{prop}}$; en caso contario, $U\geq \alpha(x_{\text{prop}}|x_{i-1})$ tome $x_i=x_{i-1}$ y $y_i=y_{i-1}$
A lo largo del algoritmo anterior es recomendable calcular el mínimo de $f$ sobre la cadena $x_0,x_1,\ldots , x_N$ construida.
# Función para propuesta con cambio por entrada
function neighbor(x₀::Array{Float64,1},i::Int64,v::Array{Float64,1})
n = length(x₀)
return (x₀ + v[i]*rand(Uniform(-1,1))*I(n)[i,:]).nzval
end
neighbor (generic function with 1 method)
# Función para temperatura en la probabilidad de aceptación
temperature(t::Real) = 1 / log(t)
temperature (generic function with 1 method)
# Algoritmo de recocimiento simulado con cambio por entrada
function recocimiento_simulado_entry(f::Function,t::Function,neigh::Function,x₀::Array{Float64,1},N::Int64,
v::Array{Float64,1}=[0.0])
n = length(x₀)
if v == [0.0]
v = ones(n)
end
x = x₀
x_v = Array{Float64,1}[]
push!(x_v,x)
x_min = x₀
y = f(x)
y_min = y
for j in 1:N
t = temperature(j)
for i in 1:n
x_prop = neigh(x,i,v)
y_prop = f(x_prop)
if y_prop <= y
x = x_prop
y = y_prop
else
p = exp( (y-y_prop)/t )
if rand(Uniform()) < p
x = x_prop
y = y_prop
end
end
if y_prop <= y
x = x_prop
push!(x_v,x)
if y_prop <= y_min
x_min = x
end
else
p = exp( (y-y_prop)/t )
if rand(Uniform()) < p
x = x_prop
push!(x_v,x)
else
push!(x_v,x)
end
end
end
if y < y_min
x_min = x
y_min = y
end
push!(x_v,x)
end
return x_v[end], x_min
end
recocimiento_simulado_entry (generic function with 2 methods)
# Función para propuesta
function neighbor(z)
[rand(Uniform(z[1] - 1, z[1] + 1)), rand(Uniform(z[2] - 1, z[2] + 1))]
end
neighbor (generic function with 2 methods)
# Algoritmo recocimiento simulado
function recocimiento_simulado(f::Function,t::Function,neigh::Function,x₀::Array{Float64,1},N::Int64)
x = x₀
x_v = Array{Float64,1}[]
push!(x_v,x)
x_min = x₀
y = f(x)
y_min = y
for j in 1:N
T = t(j)
x_prop = neigh(x)
y_prop = f(x_prop)
if y_prop <= y
x = x_prop
y = y_prop
else
p = exp( (y-y_prop)/T )
if rand(Uniform()) < p
x = x_prop
y = y_prop
end
end
if y < y_min
x_min = x
y_min = y
end
push!(x_v,x)
end
return x_min, x_v
end
recocimiento_simulado (generic function with 1 method)
# Función de Rosenbrock con un mínimo en (1,1)
function rosenbrock(x::Array{Float64,1})
(1 - x[1])^2 + 100(x[2] - x[1]^2)^2
end
rosenbrock (generic function with 1 method)
# Implementación reocimiento simulado para Rosenbrock
recocimiento_simulado_entry(rosenbrock,temperature,neighbor,[0.1,0.1],10000)[1]
2-element Array{Float64,1}:
0.9867330094947029
0.9629929384896192
gif(anim_simul_anneal, "anim_simul_anneal.gif", fps = 50)
┌ Info: Saved animation to │ fn = /home/workstation/Documents/Simulacion/anim_simul_anneal.gif └ @ Plots /home/workstation/.julia/packages/Plots/lzHOt/src/animation.jl:104
gif(anim_simul_anneal, "anim_simul_anneal.gif", fps = 50)
┌ Info: Saved animation to │ fn = /home/workstation/Documents/Simulacion/anim_simul_anneal.gif └ @ Plots /home/workstation/.julia/packages/Plots/lzHOt/src/animation.jl:104
Los algoritmos genéticos son algoritmos de optimización vía busqueda aleatoria que hacen uso de una analogía con la manera en que los genes evolucionan. Para esto se consideran tres componentes:
1) Entrecruzamiento: Sea $\mathcal{C}$ el conjunto sobre el que deseamos optimizar. Dados dos puntos $c_1$, $c_2\in \mathcal{C}$, que son interpretados como cromosomas se tiene una función de entrecruzamiento, digamos $cross:\mathcal{C}\times \mathcal{C}\rightarrow \mathcal{C}$, que produce un nuevo punto que interpretamos como el "decendiente" de $c_1$ con $c_2$.
2) Mutación: Dado $c\in\mathcal{C}$ decimos que una mutación de $c$ está dada por $mut(c)$ con $ mut :\mathcal{C}\rightarrow \mathcal{C}$.
3) Aptitud: Dado $c\in\mathcal{C}$ decimos que tiene aptitud (fitness) $f(c)$ con $ f :\mathcal{C}\rightarrow \mathbb{R}$.
4) Selección: Dada una población $\pmb{c}\in\mathcal{C}^n$ podemos se seleccionar una subpoblación $s(\pmb{c})$ con $s:\mathcal{C}^n \rightarrow \mathcal{C}^m$, $m<n$, en la que usualmente se eligen con alta probabilidad los elementos de la población inicial más aptos.
# Ejemplo del vendedor viajero
function ciudades_sim( n, map_lim)
c = Dict{String,Real}[]
for c_run in 1:n
push!(c, Dict(
"id" => c_run,
"x" => round(rand() * map_lim),
"y" => round(rand() * map_lim) ) )
end
return c
end
# calling generate_cities method to create cities
ciud_dict = ciudades_sim(5,500)
5-element Array{Dict{String,Real},1}:
Dict("id" => 1,"x" => 228.0,"y" => 131.0)
Dict("id" => 2,"x" => 94.0,"y" => 415.0)
Dict("id" => 3,"x" => 127.0,"y" => 38.0)
Dict("id" => 4,"x" => 278.0,"y" => 273.0)
Dict("id" => 5,"x" => 458.0,"y" => 343.0)
function crossover( t1::Array{Int64,1}, t2::Array{Int64,1}, cross_i::Int64)
tnew_1 = t1[1:cross_i]
for c in tnew_1
if c in t2
c_loc = findfirst( i -> i == c, t2)
splice!( t2, c_loc)
end
end
return vcat( tnew_1, t2)
end
crossover (generic function with 2 methods)
crossover( [1,2,3,4,5], [5,3,1,2,4], 2 )
5-element Array{Int64,1}:
1
2
5
3
4
function mutate(t::Array{Int64,1})
mut_1 = rand(1:length(t))
mut_2 = rand(1:length(t))
t[mut_1], t[mut_2] = t[mut_2], t[mut_1]
return t
end
mutate (generic function with 2 methods)
mutate( [5,3,1,2,4] )
5-element Array{Int64,1}:
5
3
1
2
4
function fitness(c::Array{Int64,1},D::Array{Dict{String,Real},1})
fit = 0
c = vcat(1, c, 1)
for loc in 1:length(c) - 1
p1 = [ D[c[loc]]["x"], D[c[loc]]["y"] ]
p2 = [ D[c[loc+1]]["x"], D[c[loc+1]]["y"] ]
fit += euclidean(p1,p2)
end
return fit
end
fitness (generic function with 2 methods)
fitness( [5,3,1,2,4], ciud_dict)
1597.1842840869992
function inic_pop(n_pop::Int64,l_pop::Int64,D::Array{Dict{String,Real},1})
c_vec = Dict{String,Any}[]
for _ in 1:n_pop
c = [1:l_pop;]
for i in 1:length(c)
u = rand(1:l_pop)
c[i], c[u] = c[u], c[i]
end
push!(c_vec, Dict( "id" => c, "fit" => fitness(c,D) ) )
end
return c_vec
end
inic_pop (generic function with 1 method)
inic_c_pop = inic_pop(3,5,ciud_dict)
3-element Array{Dict{String,Any},1}:
Dict("id" => [5, 4, 2, 1, 3],"fit" => 1326.9704458779634)
Dict("id" => [2, 3, 4, 5, 1],"fit" => 1477.7306971101777)
Dict("id" => [4, 5, 2, 3, 1],"fit" => 1230.4671534526028)
function evolve( inic_c::Array{Dict{String,Any},1}, n_ev::Int64, cross_i::Int64, D::Array{Dict{String,Real},1})
for i in 1:n_ev
n_pop = length(inic_c)
c_dict = inic_c[i]
u1 = rand(1:n_pop)
u2 = rand(1:n_pop)
t1 = copy( inic_c[u1]["id"])
t2 = copy( inic_c[u2]["id"])
new = crossover( t1, t2, cross_i)
new = mutate(new)
push!(inic_c, Dict( "id" => new, "fit" => fitness(new,D) ) )
end
return sort!(inic_c, by=x -> x["fit"], rev=false)
end
evolve (generic function with 1 method)
evolve(inic_c_pop, 50, 2, ciud_dict)
57-element Array{Dict{String,Any},1}:
Dict("id" => [4, 5, 2, 3, 1],"fit" => 1230.4671534526028)
Dict("id" => [4, 5, 2, 3, 1],"fit" => 1230.4671534526028)
Dict("id" => [4, 5, 2, 3, 1],"fit" => 1230.4671534526028)
Dict("id" => [1, 4, 5, 2, 3],"fit" => 1230.4671534526028)
Dict("id" => [3, 2, 4, 5, 1],"fit" => 1254.0912080872085)
Dict("id" => [1, 5, 4, 2, 3],"fit" => 1254.0912080872085)
Dict("id" => [3, 2, 4, 5, 1],"fit" => 1254.0912080872085)
Dict("id" => [5, 4, 2, 3, 1],"fit" => 1254.0912080872085)
Dict("id" => [3, 2, 4, 5, 1],"fit" => 1254.0912080872085)
Dict("id" => [1, 3, 2, 4, 5],"fit" => 1254.0912080872085)
Dict("id" => [5, 4, 2, 3, 1],"fit" => 1254.0912080872085)
Dict("id" => [3, 4, 5, 2, 1],"fit" => 1294.8367593973094)
Dict("id" => [2, 5, 4, 3, 1],"fit" => 1294.8367593973094)
⋮
Dict("id" => [3, 4, 5, 1, 2],"fit" => 1550.6099349009326)
Dict("id" => [3, 4, 5, 1, 2],"fit" => 1550.6099349009326)
Dict("id" => [3, 4, 5, 1, 2],"fit" => 1550.6099349009326)
Dict("id" => [4, 3, 5, 2, 1],"fit" => 1565.0505976063453)
Dict("id" => [4, 3, 1, 2, 5],"fit" => 1565.0506101258716)
Dict("id" => [2, 4, 3, 5, 1],"fit" => 1588.674652240951)
Dict("id" => [5, 1, 2, 4, 3],"fit" => 1588.6746647604773)
Dict("id" => [5, 1, 3, 4, 2],"fit" => 1588.6746647604773)
Dict("id" => [3, 4, 2, 1, 5],"fit" => 1588.6746647604775)
Dict("id" => [5, 3, 1, 2, 4],"fit" => 1597.1842840869992)
Dict("id" => [2, 1, 4, 3, 5],"fit" => 1820.8237731099682)
Dict("id" => [5, 3, 4, 1, 2],"fit" => 1820.8237731099684)
fit_run = Inf
c_run = Int64[]
for a in 1:5
for b in deleteat!( [1:5;], a)
for c in deleteat!( [1:5;], sort([a,b]))
for d in deleteat!( [1:5;], sort([a,b,c]))
for e in deleteat!( [1:5;], sort([a,b,c,d]))
fit_loc = fitness( [a,b,c,d,e], ciud_dict )
if fit_loc < fit_run
fit_run = fit_loc
c_run = [a,b,c,d,e]
end
end
end
end
end
end
c_run
5-element Array{Int64,1}:
1
3
2
5
4
fit_run
1230.4671534526028